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This paper proposes a novel framework for delay-tolerant particle filtering that is computationally 
efficient and has limited memory requirements. Within this framework the informativeness of a delayed 
(out-of-sequence) measurement (OOSM) is estimated using a lightweight procedure and uninformative 
' measurements are immediately discarded. The framework requires the identification of a threshold that 

separates informative from uninformative; this threshold selection task is formulated as a constrained 
optimization problem, where the goal is to minimize tracking error whilst controlling the computational 
requirements. We develop an algorithm that provides an approximate solution for the optimization 
' problem. Simulation experiments provide an example where the proposed framework processes less 

than 40% of all OOSMs with only a small reduction in tracking accuracy. 
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^ ; I. Introduction 

Tracking is frequently performed using multiple sensor platforms, with measurements being relayed to a 
central fusion site over a wireless network. This can lead to some measurements being delayed through 
packet losses or processing delays. The fusion centre is then faced with out-of-sequence measurements 
(OOSMs). For some highly non-linear tracking tasks, the particle filter significantly outperforms the 
Extended or Unscented Kalman Filter (EKFs/UKFs). Incorporating delayed measurements into a particle 
filter in an efficient manner can be a challenging task. The goal is to retain tracking accuracy while 
minimizing storage and computational requirements. 

In this paper, we propose a novel framework for delay-tolerant particle filtering that is computationally 
efficient and has limited memory requirements. To derive the framework we formulate a constrained 
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optimization problem of selectively processing only the most informative OOSMs (those that provide 
the most reduction in tracking error), where the constraint specifies a maximum allowable average 
computational expenditure. We develop an algorithm that addresses an approximation of this optimization 
problem. The method combines a Gaussian approximation of the current particle filter distribution and 
a linearization of the dynamics (similar to the EKF) to derive a procedure for rapidly predicting the 
anticipated mean squared error reduction associated with processing each OOSM. We then derive a 
threshold for selecting the "best" OOSMs while respecting the average processing cost constraint. Any 
measurements which are deemed insufficiently informative are thus immediately discarded. 

We report simulation results for an example tracking scenario where the proposed algorithm processes 
only 40% of all delayed measurements. The algorithm achieves an accuracy that is almost equivalent to 
that achievable by re -running the particle filter each time a delayed measurement is received, but reduces 
the computational cost by a factor of almost two. 

A. Related Work 

There has been substantial work on the efficient incorporation of out-of-sequence measurements OOSMs 
in Kalman filters CD-IE]]. Fewer techniques have been proposed for processing delayed measurements 
using particle filters. In ifTOl . Orton et al. propose an approach that stores sets of particles for the last £ 
time steps, where t is the predetermined maximum delay. The algorithm samples new particles at the time 
step of the delayed measurement and uses these to update the current particle weights. This method was 
improved with a Markov chain Monte Carlo (MCMC) smoothing step to mitigate the potential problem 
of degeneracy in [11]. When a large number of particles is needed for accurate tracking, the algorithm 
has an excessive storage requirement. 

Mallick et al. propose an approximate OOSM particle filter based on retrodiction in lfT2l . When 
the filter receives an OOSM, it retrodicts (predicts backwards) the particles to the time step of the 
delayed measurement and uses these particles to update the current weights. The algorithm in ifTBl also 
uses retrodiction, but employs the Gaussian particle filter of Ifl4"1 . Retrodiction requires a backwards 
information filter, i.e. a filter that runs backwards in time. Constructing such a filter is possible for 
linear state dynamics, and these are the systems that are studied in lfT2l . 1131 . Recent advances in 
particle smoothing lfl31 - lfT7l can be adopted to extend the applicability of these techniques to non-linear 
systems. However, running the backwards information filter remains a computationally intensive exercise, 
equivalent to re-running the particle filter from the time of the delayed measurement. 

In |fT31 , Orguner et al. develop strategies to reduce both the memory requirements and computational 
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complexity of OOSM particle filters. They propose a "storage efficient particle filter" that only stores 
statistics (single mean and covariance) of the particle set, rather than the particles themselves, at previous 
time steps. Auxiliary fixed point smoothers are then employed to determine the likelihood of the delayed 
measurement conditioned on each particle in the current set, and this likelihood is used to update the 
weight of each particle. The algorithm can only adjust particle weights, not change particle locations; 
this can lead to a particle degeneracy problem if an OOSM is highly informative and should induce a 
significant change in the filtering distribution. Orguner et al. propose a heuristic of ignoring OOSMs that 
lead to filter degeneracy, but this is not satisfactory, since the highly informative OOSMs are often the 
most important to process. 

The algorithm we propose in this paper involves selective processing of OOSMs. This was first 
discussed by Orton and Marrs in |[T0l ; they advocated a heuristic approach of discarding all measurements 
that are delayed beyond a constant time, with the constant to be determined through experiment. More 
recently, selective OOSM processing has been considered by Tasoulis et al. in |[T8l and in our previous 
work |fl9ll . Tasoulis et al. proposed a number of heuristic metrics to estimate the utility of delayed measure- 
ments and develop threshold-based tests to discard measurements of low utility. They incorporate these 
tests into three Kalman filtering algorithms that are designed to process delayed measurements. In fl9ll 
we proposed a threshold based procedure to discard uninformative delayed measurements, calculating 
their informativeness using mutual information and Kullback-Leibler distance metrics. We applied our 
approach in the general non-linear setting, using a combination of the storage-efficient particle filter 
proposed in |[T5l and a re -run particle filter. 

The approach proposed by Tasoulis et al. is developed for the Kalman Filter and it is difficult to extend 
to more general filtering problems with non-linearities. The proposed utility metrics are heuristic and do 
not truly capture the potential that each delayed measurement has to improve the tracking performance. 
The latter issue is also a failing of our own work in |fT9ll ; although mutual information and Kullback- 
Leibler distance metrics measure the potential for information gain, they do not directly assess the potential 
reduction in estimation error. Perhaps most importantly, neither |fl"8l nor (I9\i identifies a procedure 
for threshold selection, despite the fact that the choice of this threshold can have a major impact on 
performance and the appropriate value is a highly application-sensitive quantity. 

B. Paper Organization 

The rest of the paper is organized as follows. Section HTl provides a formal problem statement. Section UlTl 
describes memory efficient OOSM particle filters. Section [TV] presents the proposed novel framework 
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for selecting informative OOSMs. In Section [V] we explore the approximations made in the derivation 
of the framework and present a theorem identifying asymptotic conditions under which one of the key 
approximations becomes exact. Section [VT] presents a concrete OOSM particle filtering algorithm based 
on the selection framework and Section IVIII describes simulation experiments for an example tracking 
scenario. We make concluding remarks in Section IVIII I 

II. Problem Statement 

We now provide a formal statement of the OOSM filtering problem that we address and formulate the 
optimization task. We consider the general discrete-time Markov state-space model with state dynamics 
and measurement models both defined by non-linear maps. The innovation and observation noises are 
modelled as additive Gaussian. At each timestep k, there is an active set of distributed sensors, V k , that 
make measurements and K = sup fc>1 |Vfe| is the maximal number of active sensors. These measurements 
are relayed to the fusion centre. A subset of them Sk experience minimal delay and can be processed 
at time k. Other measurements are delayed and only become available for processing at later timesteps. 
Measurements delayed by more than I timesteps are discarded. 
The system is described by the following state-space model: 

X k =f k (X k _ 1 ) + ti k (1) 

Y k s = h s k (X k ) + C k (V S GV fc ) (2) 

y k = {Y k Sk : S k C V k } (3) 

^k — X^k-i i x k-l+l '•••» 1 fc-l J W 

Here {X k } denotes the state sequence, which is a Markov diffusion process with initial distribution 
Xq ~ p{xq), and {Y k s } denotes the measurement sequence at the s-th sensor, with Y k Sk = {Y k s : s £ S k }. 
■dk is the innovation noise with Gaussian distribution M(0,V k ), and Q is the measurement noise with 
Gaussian distribution A/"(0, Q|). The functions f k : R d -> R d and hf. : R d -> R m ° are the state transition 
and measurement maps. 34 denotes the set of non-delayed measurements received at time k. Z k denotes 
the set of OOSMs received at time k. The set S Tjk is the subset of active sensors at time r whose 
measurements are received at time step k {S k ^k = Sk)> Yr T,k is the set of measurements made at time r 
that arrive at the fusion centre at time k. 

1) OOSM Filtering: Let Wi : j k denote the set of measurements generated in the interval available 
at the fusion centre by time k. This includes all the non-delayed measurements y^j = \J J m=i y m and 
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OOSMs Z i: j jk = {Z s Tm G Z m : r G [i,j],s G S Tim ,m e [i + l,k]}, where Z^ m is the OOSM that was 
acquired at time r by the sensor s G <S T ,m C V r and was received at the fusion centre at time m. Let 
YVi-.j,k = Wiy^ \ -Hfc = {Wi:j,k-l,yk}, i- e - the set of all measurements available at time k except those 
in Z k . Lastly, note that Z T>k = Z T ., T ^ k , W T:k = W T:T , k and W Tjfc = W T:Ttk . 

The sequential OOSM filtering task involves calculating an estimate X k of the current state, given 
all available measurements at time k, Wi :k . k . In this work, we form the estimate by calculating an 
approximate expectation of the state by sequentially calculating a particle representation of the posterior 
distribution. 

2 ) Selective Processing for Computational Constraints: In this paper we are interested in reducing 
computational requirements by processing only the informative OOSMs. We formulate this problem as 
an optimization problem that involves minimizing the mean-squared error (with respect to an L2 norm) 
subject to satisfying a constraint (C ave ) on the expected computation at each time step. 

Let b s rk G {0,1} be the indicator of OOSM Z s rk arrival and denote by p s rk the expected value of 
b s T k , conditioned on all the measurements received prior to time k. Denote by k the computational 
cost associated with processing the OOSM Z s rk . Let d s rk G {0, 1} be our decision to process or reject 
measurement Z s r k and T> k = [j^Z k _ i Uses k fe) ^ e the current set of all possible decisions. Decisions 
must be made sequentially, prior to the arrivals of the OOSMs at time k due to the real-time nature of 
the tracking task. The goal is to ensure that the computational constraint is obeyed on average at each 
time step, i.e. in expectation with respect to all possible arrivals of OOSMs. 

We thus address the following optimization task for each k over the tracking period: 



Previously proposed OOSM particle filters primarily differ in how they incorporate the OOSMs from the 
set Z k . The simplest approach is to discard them, but this often results in poor tracking performance. 
Another obvious approach is to restart the filter at the time step immediately prior to the time step 
associated with the earliest OOSM in Z k and re-run to the current time step k. This requires that we 
record all the particles, weights and the measurements for the maximal delay window. We call this 
approach the "OOSM re-run particle filter" and consider it to be an accuracy benchmark. This method 
has two unattractive qualities: the storage requirements can be immense and the computation cost is high. 




K k-l 



SUbjeCt tO J2 J2 d r,kPr,kCr, k < KC ; 



avc 



(5) 



fc=i T =k-i sev T , fc 



III. OOSM Particle Filters 
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As discussed in Section ll-Al several methods have been proposed to alleviate these costs. In this section, 
we provide a brief review of the storage efficient particle filter of |[T5l and describe a relatively obvious 
alternative algorithm that we introduced in |fl9l . In both algorithms, the memory requirements are reduced 
by storing statistics of the particle sets from past time steps instead of the particles themselves. The past 
particle distributions are approximated by Gaussian approximations. The stored information is then the 
mean and covariance matrix of particles at each time step from k — I — 1 to k. Denote, respectively, by 
u>k the sets of the values and weights of particles at time k, and let \x k , H k denote their mean and 
covariance. The stored information is then 

&>k = {Hk-e:k,R-k-e-.k,yVk-e-l:k,k}, (6) 

Here Hk-£-i:k an d Hk-i-v.k denote, respectively, the means and covariances of the particle sets for 

time-steps ranging from k — t — 1 to k. 

A generic storage efficient OOSM particle filtering algorithm is summarized in Algorithm Q] If there 

are no OOSMs at time k, we write Z k = 0. 
Algorithm 1: Generic OOSM Particle Filter 

1 At time k 

Input: Z k , w fe _iJ4 

2 (£fc, Wfe) <- ParticleFilter (34, £ fc _ b ; 

3 O fc , R fc ) <- SaveGauss (£ k , Vk) 

4 if Z k ^ then 

5 (£ fc , u> k , ri fc ) <r- ProcessOOSM(Z fc , £ fc , u k , ft k ) ; 

In this algorithm, the function ParticleFilter can be any standard particle filtering method. If 
y k = 0, ParticleFilter only propagates the particles and skips the measurement processing step. 
The function SaveGauss calculates the maximum likelihood estimates of the mean and covariance 
given the weighted sample set £ k , u k and stores these in Q k : 

N 

M* = E"M (7) 

1=1 

N 

1=1 

The function ProcessOOSM specifies how OOSMs are processed and varies depending on the specific 
algorithm. 
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A. Gaussian Approximation Re-run Particle Filter (OOSM-GARP) 

A simple modification of the re-run particle filter involves storing only Gaussian approximations of the 

particle distributions at previous timesteps. When a batch of OOSMs arrives, the particle filter is re-run 

from the time step preceding the earliest OOSM. Since the particle set from that time step is unavailable, 

particles are generated from the stored approximation. 

When OOSM-GARP receives Z k at time k, it returns to the time step r k — 1 (let r k denote the earliest 

time step of all OOSMs in Z k ). It samples particles from A/"(/x^_i, R^-i), propagates them to the 

time step r k and runs the filter as standard particle filter using all stored measurements. At each step, it 

updates the mean and co variance matrix in the stored set ft k as described in Algorithm [2] 
Algorithm 2: ProcessOOSM-GARP 
Input: Z k , fl k 

l T fc = min{r : y T G Z k } ; 

T 

3 41-i = W, i = i...N ■ 

4 for j = r k ,...,k do 

5 uj) <r- ParticleFilter (W jtk , Uj-i) ; 

6 (/Xj, Rj) <— SaveGauss cjj) ; 

7 endfor 

In many tracking tasks, the Gaussian provides a reasonable approximation to the particle distributions. 
In OOSM-GARP, the Gaussian is only used to re-start the particle filter (to draw initial samples), so the 
impact of approximation errors on filtering performance is relatively small. OOSM-GARP thus performs 
almost as well as the basic re-run particle filter but requires much less memory. However, OOSM-GARP 
is relatively computationally complex since it reprocesses all the particles for k — r k + 1 steps. Note that 
the cost to process OOSMs corresponding to a single time step, Zj jk , is approximately equal to that of 
processing the whole batch of OOSMs Zj :k -i tk since we have to execute the particle filter from time 
j to time k in both cases. This cost is proportional to the total computational complexity of functions 
ParticleFilter and SaveGauss multiplied by a factor of k — j. 

B. Storage Efficient Particle Filter with EKS (SEPF-EKS) 

We now provide a brief review of the storage efficient OOSM particle filter from |[T5l . Orguner et al. 
described three versions of the filter, which differed according to the auxiliary fixed-point smoother 
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they employed. We focus on the filter that employs Extended Kalman Smoother, since it is the least 
computationally demanding but has comparable tracking performance. 
The SEPF is based on the following weight-update equation: 

K P (Zr tk \$,W 1:ktk )J&. (9) 

(i) (i) 

Here u k L and u k denote the weights before and after processing Z T k . The SEPF estimates this like- 
lihood expression in two stages. First it approximates p{x T \£ k \Wi :kjk ) by applying an augmented-state 
extended Kalman smoother ll20l . treating the current particle £^ as a measurement. The SEPF then em- 
ploys an EKF approximation of p(Z Tk \x T ) to construct an estimate of the likelihood p(Z Tk \£ k , W\ :k k ). 

Although the original algorithm was designed to treat individual OOSMs, it can be easily extended 
to treat batches of OOSMs by running a separate update for each time-step. This extended algorithm is 

presented as Algorithm [3] 
Algorithm 3: ProcessSEPF-EKS 
Input: Z k , Q k , u} kjf , £ k 

1 for Z T)k G Z k do 

2 Compute approximation p{Z Tjk \£ k , Wi: ktk ) for all i; 

3 <r- u^p(Z Tjk \$\w 1:kjk ) Vi ; 

4 endfor 

SEPF-EKS achieves significant computational savings because the filtering operations for step 2 are 
common to all N particles except for a single time-step. This means that the effective computational cost 
is equivalent to running one time step of a particle filter, and is therefore usually less than that of the 
OOSM-GARP filter. The advantage diminishes when it is common for OOSMs to arrive in batches with 
different delays because of the seemingly unavoidable loop in the algorithm. 

IV. Selective OOSM Processing 

The computational cost of processing an OOSM is relatively high and frequently it is wasted effort, 
resulting in minimal change to the filtering distribution or the tracking accuracy. In this section we design 
a procedure for addressing the optimization problem posed in Section UU that of minimizing the mean 
squared error while controlling the computational effort. 

The optimization problem is challenging and generating an exact solution would be more costly than 
simply processing all OOSMs with a re -run particle filter. We therefore strive to approximate the problem 
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so that we can develop an efficient procedure for selecting the informative OOSMs. The complexity of 
this procedure must not depend on the number of particles in the filter. 

Our method employs a Gaussian approximation of the joint distribution of the current state and the 
current set of OOSMs. We derive this approximation using an EKF-type linear approximation of the 
general state-space model. Second, we model the OOSMs from different sensors or different times 
as approximately unconditionally independent. This second approximation allows us to disentangle the 
effects of processing different OOSMs on the filtering error. In section [V] we study asymptotic conditions 
under which this assumption holds exactly. This provides a solid theoretical justification for our choice of 
this simplifying approximation and we consider that it is sufficiently accurate in practice for our purpose 
of selecting the informative OOSMs. It is important to stress that these approximations are only used for 
the purpose of selecting the measurements to process; they are not employed within the filter itself. 

A. Tracking MSE Under Gaussian Approximation 

We employ the well known EKF-type linear approximation of the general state-space model: 

X k = fkinx^) + F fc (X fc _! - nx^) + $ k (10) 

Yk = K^x k ) + n s k (x k - » Xk ) + ci, s e v k (ii) 

Here F k and H| are linearizations (through Taylor expansion at ftx t l and /j,x k , respectively) of the 
non-linear dynamic and measurement maps. 

Let Pfc and \x k be the covariance matrix and the mean of the Gaussian approximation of the joint 
probability distribution of the current state and the current set of OOSMs conditioned on all available 
measurements. The covariance matrix and the mean have the following structure: 

R x fc x fc |w 1:fc , fc ^x k z k \w 1:k , k I = I ^x k \m :k , k I q 2 ) 
^z k x k \m.. k , k R 2 fc 2 fc |w 1:fc , fc / \^z k \m..KkJ 

where R„ „ is the current state covariance, R„ „ ,.~ and R„ „ = R T ~ 

AjiAfel KVl : fc,fc Ajfc-Sjfel Wl: klk -SfcAfc | Wl- k , k X k Z k \V\)\- kk 

is the state-measurement cross-covariance and R, _ ,,-77 is the measurement set covariance. Note 
that the means and covariances are conditioned on Wx : k,k which includes the current set of undelayed 
measurements y k as well as all the OOSMs and undelayed measurements that have been incorporated 
up to time k. In the following discussion, we will often skip this conditioning to avoid unnecessarily 
complicated notation, but this conditioning is implied unless explicitly stated otherwise. 

The optimal MMSE estimator Xu of the state is known to be the conditional mean u v ,cr, ~ , 

Afc I Wl: k . k ,<£ k 
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which in the case of our Gaussian approximation is simply: 

x k = Vx k + R-x k z k R-zlz k (Zk - Hz k )- (13) 

Let B k = Ur=fc-^Uses k iK k\ b e ^ e set °^ ran dom variables that indicate OOSM arrivals at time k. 
This set defines the structure of the set Z k along with the associated mean nz k and (cross-)covariance 
terms H Xk z k and Hz k z k - By the law of total variance the variance of the the estimator can be expressed 
as: 

var(X fc - X k ) =E{var(X fc - X k \B k )} + var(E{X fc - X k \B k }). (14) 

Since, according to our linearization, E{X k \B k } = A(/ix t _J and E{X k \B k } = [i Xk = /fcO^.J we 
have for any realization of B k : E{X k — X k \B k } = 0. Thus the variance of the MMSE estimator is equal 
to the expectation of its variance conditioned on the realization of indicators B k : 

var(X fc - X k ) = E{var(X fe - X k \B k )}. (15) 

For a specific realization of indicators B k this variance is defined by the components of the joint covariance 
matrix (recall that Z k is a function of B k ): 

var(X fc - X k \B k ) = R XkXk - Rx k z k ^l z Rz k x k (16) 

The mean squared error of estimating the state X k conditioned on the OOSM set Z k (as well as all the 
previous measurements) is thus given by 

trvar(X fc - X k ) = E{trvar(X fc - X k \B k )} (17) 
= tr K XkXk - E{tr R M R^R M , } (18) 

Under the assumption that the measurements made by different sensors (or the same sensor at different 
times) are approximately unconditionally independent, Hz k z k is approximately block-diagonal. This 
implies that we can approximate the above expression as follows: 

trvar(X fe - X k ) « tr K XkXk -E j d *,kK,k tr R^y-R^Ry-^ } . (19) 

Here the expectation is taken with respect to the measurement arrival indicators b s rk , RyY s = var(Y!j?) 
is measurement covariance and Ry s x k = cov(Y r s ,X k ) is the state-measurement cross covariance. If we 
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denote 

R t,k = tr R X fc Y/ Ryjy^s ^Y*X k > (20) 

the factor that we will refer to as the measurement utility then the expression for the MSE can be further 
simplified: 

fe-i 

trvar(X fc -l fc ) «trR Xfc x fc - ^ £ <Xfc^,fc- (21) 

T= k-ls&V T 

where p s rk = E j 6* k j is the probability that the measurement acquired by sensor s at time r arrives 
at time k (conditioned on the measurement arrivals up to time k). The above expression is a natural 
objective function to be minimized to assure the best tracking quality. The minimal value of the objective 
is reached when all measurements are processed (d s rk = l,Vr, s) since R^ k >0. 

B. One-step Constrained Minimization of Approximate MSE 

Given the discussion above and the identified approximations, the constrained optimization problem 
posed in Section HT1 can be formulated as follows: 

min tivax(X k — X k ) 

£> fc e{0,l} E - ' v - 

fc-1 

subject to £ d r,kPr,kCr,k < C ave (22) 

T=k-£s€V Tik 

The unconstrained objective to be minimized can be formulated using Lagrange relaxation with Lagrange 
multiplier 7^: 

J(V k ) = tr var(X fc - X k ) + 7 fc ( ]T £ < )fc p^ |fc C' fc - C ave J 

Yr=fc-^seV T / 

fc-1 

= trR XfcXfc - £ E <*Pr I *(^,fc-7fcC? ifc )-7*C'«vB- (23) 
r=fc-£ sev T 

For a fixed 7^ the optimal solution can be found by optimizing each d s rk independently since the 
contribution of each term under the sum corresponding to a particular d s T k is independent of all other 
variables to be optimized. It is clear that setting d s rk = 1 whenever R s rk — j k C^ k > and d s T k = 
whenever R s rk — jkC^ k < produces the smallest value of the objective function for a given 7^. 
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Substituting this solution into the constraint we obtain 

fc-i 

2 2 l {R s r:k ~lkC^ k }P S r,kCr,k < ^ave, (24) 

T=fc-^SgV T 

where is the indicator function. If we denote R s rk = R s rk /C^ k , the measurement utility diminished 
by the processing cost incurred, the above is equivalent to 

Pr,k^T,k ^ Cave- (25) 

The optimal value of 7*. is thus the smallest value for which (T25T ) holds. A simple practical algorithm 
can be devised to identify this value of 7^. The algorithm, summarized in Algorithm HJ assumes that we 
can evaluate p s T k , which is usually possible given sufficient knowledge about the measurement apparatus 
and the network delay profile. 

Algorithm 4: Threshold selection algorithm 
Input: {R s T k }, IX J, {C s Tk } of cardinality T = E^iL \ V r\ and C ave ; 

1 Order set {R s Tk \ by decreasing value, output ordered sequence {i?°}^ =1 ; 

2 Construct sequences {p°}^ =1 , {C°}^ =1 using mapping (r, s) 1-4 n used for the previous set ; 

3 Construct sequence { , f°}^ =1 with elements = Ej=iPjC° 

4 Identify n* = arg max n : ^° < C a ^ e 

Output: 7 fc = ; 

We can now describe the operation of the proposed OOSM selection algorithm. At every filtering step 
the selection algorithm first calculates the measurement utilities diminished by the processing cost, R s rk , 
along with probabilities of arrival for all possible OOSMs, p s rk . It then identifies a threshold 7^. such 
that the expected processing cost does not exceed C ave (step 4 in Algorithm @]). The final step of the 
algorithm is to select arriving OOSMs with utility R s r k surpassing the calculated threshold. 

To execute the proposed algorithm we need expressions for the (cross-) covariance matrices Hx k Y T s and 
Ry s y». These matrices can be calculated online using the extended Kalman smoother (EKS) algorithm. 
We employ the Rauch-Tung-Striebel (RTS) backward recursion realization ETTl . We apply the RTS 
recursion starting from the Gaussian approximation of the posterior at the current time k and moving 
backwards in time until time step k — I. As a result, we obtain a sequence of smoother means [i x ^ 
and covariance matrices Hx T x T for k — £ < t < k. 

At time k we have the set of measurements VVi : fc,fc, so the linearizations (ITOb can be made more general 
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(and, hopefully, accurate) with the use of the EKS statistics fj, x ^ < k — 1: 

y T s = K(i*x T \m.. h . h ) + H -( x - " ^m-.J + £, * e v T . (27) 

Here the Jacobians F T and are evaluated at the points defined by the respective EKS means. With 
the use of the above linearization, calculation of the required approximate covariance matrices becomes 
straightforward. Noting that E{Y r s } = h s T {yix T ), and observing the independence of Q and X T — [ix T , 
we can derive 

Ry-y- = H s t R Xt xH s t T + R C?C f (28) 

Note that Rx T x T is the covariance of the extended Kalman smoother. 

Next, we calculate the cross-covariance Rx h v>- Since E{X r } = f T (^x T ^ 1 ), we have for any r < k: 

X k -E{X k } = F k (X k _ 1 -pL Xk ^) + # k (29) 
= F ifc (F fc _ 1 (X jt _ 2 - MXfe _ 2 ) + ^_ x ) + $ k (30) 

= F fc F fc _i(F fc _ 2 (X fc _ 3 - ^x fc _ 3 ) + #k-2) + Ffc^fc-i + $k (31) 

k 

= F k>T (X T - nxj + F kj#j (32) 

j=T+l 

where we have introduced the notation F k r = Yl J=T+1 Fj and F k>k = I. We can thus evaluate the 
cross-covariance using the expression: 

Rx k Y T ° = F k!T Rx T x T FL s T T . (33) 

V. Asymptotic Optimality of the Proposed Algorithm 

In this section we will consider the conditions under which the unconditional measurement inde- 
pendence approximation made in the previous section is expected to hold, assuming that the Gaussian 
approximation is accurate. The assumption simplifies the algorithm derivation and reduces its computa- 
tional requirements, but it leads to sub-optimality of the derived constrained MSE minimization algorithm. 
The conditions established in this section help us understand when the performance of the proposed sub- 
optimal algorithm is expected to approach that of the optimal OOSM selection algorithm, assuming that 
the Gaussian approximation and linearization are accurate. 
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The following theorem specifies that, under mild regularity assumptions, if an asymptotic condition on 
the minimal eigenvalues of the noise matrices holds, then the block-diagonal approximation employed 
to derive the OOSM selection algorithm in the previous section holds exactly. The proof is provided in 
Appendix |A] 

Theorem 1. Let ~Rz k z k , Rx t 2 fc be defined as in (1121 ) and let ~Bz k z k be the block-diagonal matrix whose 
diagonal blocks match those of ~Rz k z k (the covariances of measurements from the same sensor at the 
same time). Suppose that the following assumptions hold: 
A\. p(Rx n xJ < oo, \fk - £ < n < k 

A 2 : p{W m W m T ) 1 / 2 < oo and p{W m T W m ) 1 / 2 < oo, VA; - I < m < k and Vs G V m 
A3: pfFm^F^) 1 / 2 <oo, \fk - £ < n < m and k - £ < m < k 
Then we have for any £, K < 00, and k > 1: 

mm A min (R aa ) ->■ 00 | tr R Xk z k ~R-zl Zk ~Rz k x k - tr Rx^.B^^R^x, | 0, (34) 

where A min (-) = min, Aj(-) 

The regularity conditions imposed in Theorem Q] are mild and natural. Assumption Ai requires the 
extended Kalman smoother covariance Rx„x„ to have finite spectral radius. Thus assumption A\ basically 
requests the stability (including the numerical stability) of the EKS. Assumption A2 and .A3 require 
the spectral radia of matrices H^H^ T , H^ T H^ and F m n F^ n to be finite. If the measurement and 
transition functions, ht{-) and /&(•), are differentiable (sufficiently smooth), leading to F m and with 
finite elements, then assumptions A2 and A3 hold by the Gershgorin disc theorem |[22l . Any scenario 
when the EKS functions normally and can be implemented leads to the assumptions being satisfied. 

The asymptotics in the theorem are with respect to min s m A m ; n (R^^ ) — s> 00. The implications are 
best illustrated by way of example. If R^ s £» is scalar for all sensors at all times and R^* = C7 2s m , 
then A m in (Rr» ) = o~ 2S m . The asymptotic condition thus implies that the measurement noise variance 
approaches infinity for all sensors at all times, or, equivalently, that all measurements become utterly 
uninformative. If, on the other hand, R^ is 2 x 2 with equal component variances, for all sensors at 
all times: 



then Amin (R^» ^ ) = a 2S m (l — Thus the asymptotic specifies that measurement components are 




(35) 
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not absolutely (positively or negatively) correlated (|r£j ^ l,Vm,s), and that they have asymptotically 
large variance. 



In this section we specify an OOSM particle filter that employs the general OOSM selection framework 
presented in Section [IV] For clarity, we describe the filter in the context of a specific application scenario, 
but it can be easily adapted to different delay models and OOSM processing costs. 

We consider a situation when there are several sensors sending measurements (e.g. bearing or range) 
of the target to a common fusion centre. All sensors are assumed to have communication issues leading 
to OOSMs. An OOSM arrives at the fusion centre from a given sensor with probability p osm and delay 
d. The delay d is uniformly distributed in the interval [0,£]. The probability 1 — p osm characterizes the 
events that (i) the OOSM is dropped in the network; reaches the fusion centre at all. For example, this 
corresponds to the scenario when OOSMs delayed by more than i are automatically dropped by the 
network. 

We implement the proposed OOSM processing framework using the SEPF-EKS algorithm of |fT51 . 
In this case the OOSMs with the same time stamp arriving from different sensors can be processed 
in one sweep of SEPF-EKS algorithm (see Algorithm [3] and associated discussion). Instead of a single 
OOSM Z s rk we thus consider a set ZJ:/ consisting of {Z^_ k } and designated by the ordered index set 
X = {0, lj'^l such that Z s rk ^ zffi if and only if the element corresponding to sensor s, Z(s) = 1. 

We set the cost to process Z^ fc as C;^ J = 1, (the cost to run the SEPF-EKS algorithm on a given 

•fx) 

hypothetical realization Z^/), irrespective of the particular combination of Z s rk . We make this choice 
because the computational complexity of the SEPF-EKS algorithm is approximately the same as one 
timestep of the particle filter. The average cost constraint analogous to (l22l) is then: 



where d T I is the decision whether or not to process a given realization / and J is the set of all 
possible realizations of X. C ave can be interpreted as the average number of SEPF-EKS algorithm sweeps 
per filtering step or, in other words, the average additional overhead caused by OOSM processing. 
We have an expression for the MSE analogous to (1211) : 



VI. Selective OOSM Particle Filter 




(36) 



fc-i 



trvar(X fc - X k ) = trRx fc x, 



d r,k Pr,k u r,k ■ 



(37) 



m=k-ei£3 
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Here is calculated similarly to (l2Qb with Y m being the vector constructed from those for 



which X(s) = 1: 



R ml - tr R JC fc yi I} R yi I} x fc • ( 38 ) 



The probability £r fc J that an OOSM with a given sensor combination Z active at time r arrives at time 
k can be calculated as: 

p { 3= n n t 1 ( 39 ) 

seV T :X(s)=l jGV T :X(i)=0 

Here ■ = if the measurements from sensor s at time r have already arrived. If not, then: 

p s Tjk = Pr{A* = k - r|6' )fc _! = 0, . . . , 6* >T = 0} (40) 

where A* is the delay that the OOSM from sensor s experiences at time r. For the case of the uniform 
delay distribution and probability of successful transmission p osm , we have: 

V ^ = l + l-(k-T) (41) 

Equipped with the expressions above we can calculate RiJ, the analog of R* k , and apply a slightly 
modified version of Algorithm |4] to set the threshold 7^. This algorithm employs a similar measurement 
covariance matrix block-diagonality approximation as in the general framework described in Section |IV] 
In this case, however, the blocks are larger and consist of matrices m Y m Y {x}, rather than Ry^y^. The 
modified approximation is thus that blocks R y {z> y {j} are close to zero for all combinations of sensors 
I and J and any m 7^ n. 

As a final heuristic refinement of the algorithm we use the OOSM-GARP algorithm to process 
those OOSMs for which the SEPF-EKS algorithm performs poorly, namely the highly informative 
measurements that should induce significant shifts in the current filtering distribution. We add a test to 
check whether the effective number of samples in the particle filter drops significantly after the application 
of the SEPF-EKS processing; if this occurs, we reprocess the OOSM using the OOSM-GARP filter. This 
allows the algorithm to adjust both weights and locations of particles to account for the new information 
embedded in the OOSMs. We have observed that this step greatly improves the performance of the filter 
in difficult situations at a minimal cost. 

The OOSM particle filtering algorithm based on the above discussion is presented in Algorithm [5] 
This algorithm describes only the OOSM processing procedure corresponding to ProcessOOSM in 
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Algorithm Q] (Algorithm Q] presents the complete high level OOSM particle filter pseudocode). In Algo- 
rithm [5j as we discussed in Section [TV] we first calculate the sequence of EKS means and covariance 
matrices, which are further used to compute the Jacobians and the utilities {R T ^}. These are used in 
CalcGamma (a minor modification of Algorithm H}, which has the task of setting the current value of 
threshold 7^. This threshold is used to determine which OOSMs should be processed with the function 
ProcessSEPF-EKS summarized in Algorithm [3j The failure of this algorithm, which is expressed 
through particle degeneracy, is detected via the second threshold test (where the value of v should be 
small, e.g. 1/40). If a failure is detected, the algorithm switches to recalculate the current particle set via 

function ProcessOOSM-GARP summarized in Algorithm [2] 
Algorithm 5: Particle Filter with selective OOSM processing (ProcessOOSM) 

1 At time k 

Input: Z k , u) k , fife, C ave 

2 (VX k _ lk , Rx 4 _ t:t X k -«) <~ EKS (Q k ) ; 

3 (7 fc , {R T ,k}) <~ CalcGamma Ri t _ M x t - (:t . Cave) i 

4 EKSfailed = ; 

5 for r : Z Tjk G Z k do 

6 if R Ti k > jk then 

7 (JVpriar) ^- l/||Wfc||l ; 

8 (Wjfc) <- ProcessSEPF-EKS (Z T;fc> £ k , oj k , ft k ) ; 

9 (A^post) Vllwfclll ; 

10 if A^post < z;Ap r i or then 

11 EKSfailed = 1 ; 

12 break ; 

13 endif 

14 endif 
is endfor 

16 if EKSfailed then 

17 (cj fc , £ fc , fi fc ) ProcessOOSM-GARP (Z h fl k ) ; 



VII. Numerical Experiments 

In our simulations we consider a two-dimensional scenario with a single target that makes a clockwise 
coordinated turn of radius 500m with a constant speed 200km /h. It starts in the y-direction with initial 
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position [—500m, 500m] and is tracked for 40 seconds. 

The target motion is modeled in the filters by the nearly coordinated turn model f2~3l with unknown 
constant turn rate and cartesian velocity. The state of the target is given as Xk = {pf.,p k ',V);,v k ',u>k\ T , 
where p, v and to denote the position, velocity and turn rate respectively. The dynamic model for the 
coordinated turn model is 



(i 





sin(wfc) 


cos(cj fc )-l 


o\ 


OJfc 


Ulk 





1 


l-cos(aj fc ) 


sin(w fc ) 







Wfe 










cos(w fc ) 


- sin(a;fc) 











sin(a; fc ) 


cos(u k ) 





v° 











V 



Xk + flk+i 



where $ k +i is Gaussian process noise, ~ 7V(0,V fc+ i), V fc+ i = diag([30 2 , 30 2 , 10 2 , 10 2 , 0.1 2 ]), 

and the sampling period is 1 second. We assume that the filter initially knows little about the state of 
the target and therefore it is initialized with the state value nx a = [0, 0, 0, 0, 0] T and a large covariance 
R XoXo = diag([1000 2 , 1000 2 , 30 2 , 30 2 , 0.1 2 ]). 

There are three sensors SI, S2 and S3 sending bearing-only measurements of the target to a common 
fusion centre. The sensor locations are [Sf,Sf] = [-200,0], [Sf,Sf] = [200,0], [Sf ,Sf] = [-750,750] 
and the bearings-only measurement function is: 



hl(xk) = arctan( 



Pk 



v _ s y 



J = 1,2,3. 



(42) 



Pk ~ *j 

The measurements from the sensors are corrupted with additive independent Gaussian noises with zero 
mean and standard deviation a s = 0.05. An OOSM arrives at the fusion centre from a given sensor with 
probability p osm and delay d. The delay d is uniformly distributed in the interval [0,5]. The probability 
Posm that an OOSM reaches the fusion centre at all is set to 0.7. 



A. Benchmarked Filters 

We have implemented five different particle filters, all based on the Sampling Importance Resampling 



(SIR) filtering paradigm 1241 . The prior distribution is used as the importance functiory. The filters were 
implemented in Matlab and the code was highly optimized. 



'Although better performance could be achieved by using a more carefully-chosen importance function, this generally comes 
at the cost of some computational expense. By using the same, simple importance function for all particle filters we achieve a 
fair performance comparison. 
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Fig. 1. Tracking performance of the particle filters as a function of time using RMS error as a performance metric, (a) The 
curves show the means of 5000 Monte-Carlo trials, (b) Errorbars showing the variation of position RMS for SEPF-EKS, PF-GS 
and PF-SEL, when they use 2000 particles. The box has lines at the lower quartile, median(red line), and upper quartile values. 
Outliers (red '+') are values beyond the range of 5 times the interquartile range from the ends of the box. 
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PFall: collects all measurements from all active sensors (no OOSMs). This is an idealized filter that 
provides a performance benchmark; a real-time implementation is impossible. 

PFmis: discards all OOSMs and therefore only processes the measurements with zero delay. 
SEPF-EKS: Storage efficient particle filter using EKS smoothing as described in lfl5l (Algorithm [3). 
PF-GS: The OOSM-GARP algorithm described in Algorithm |2] 

PF-SEL: Selective OOSM processing based on the proposed framework and described in Algorithm [5] 
We use the root mean-squared (RMS) position error to compare the performances of the particle filters. 

Let and {p%i,P v ki ) denote the true and estimated target positions at time step k for the i-th of 

M Monte-Carlo runs. The RMS position error at k is calculated as 



B. Results and Discussion 

In our first experiment we fix the computational cost C ave = 0.6. We thus allow 0.6 sweeps of the 
SEPF-EKS algorithm to be performed on average per filtering step. Over an extended period of time of 
sufficiently large length L this leads to an additional OOSM processing overhead of ~ 0.6LCeks where 
Ceks is the cost of one sweep of the SEPF-EKS algorithm. If we do not apply the proposed procedure 
and process all the available OOSMs this cost in our application scenario is approximately 1.5LCeks- 
This implies that we process only approximately 40% of all measurements. 

In Fig. [T] we plot the respective RMS position performance for the tracking period of 40s for the 
algorithms with these settings. Corresponding error-bar plots of the RMS performance are shown in 
Fig. |l(b)| The actual number of individual OOSMs processed by the SEPF-EKS after application of 
the first threshold 7*. measured in our experiment is 40.04%. After the second threshold v = 1/40 the 
percentage of most informative OOSMs processed by rerunning the particle filter using OOSM-GARP 
is 1.57%. 

Fig. CD indicates that despite processing only a relatively small fraction of the OOSMs, the proposed 
algorithm performs almost as well as the much more complex OOSM-GARP algorithm (PF-GS). The 
calculation of the selection criterion has minimal overhead, so discarding the uninformative measurements 
results in significant computational savings. Thus the proposed filter is more computationally efficient 
than the SEPF-EKS filter and yet, as can be seen from Fig. \T\ it has better RMS performance. Fig. |l(b)| 
indicates that the performance of SEPF-EKS is not as stable as that of PF-GS and PF-SEL. In the 




(43) 
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proposed algorithm the increased robustness and performance stability is achieved by using the second 
threshold to detect situations when reweighting particles induces sample degeneracy problems. 

In our next experiment we study the computational complexity versus accuracy trade-off for the pro- 
posed algorithm. We illustrate this by varying the computational complexity of the proposed algorithm by 
adjusting C ave and plotting the RMS error vs. computational load measured in MATLAB. We use the fol- 
lowing values to control the OOSM processing overhead: C avc = {0, 0.1, 0.2, 0.4, 0.5, 0.6, 0.8, 1, 1.3, 2}, 
72 = 1/40. These results are reported in Fig. [2] In this figure we show the relationship between complexity 
and performance for the proposed algorithm with ten values of C ave and results of 10 simulations for 
other algorithms. Each simulation involves 1000 Monte Carlo runs. We compare the performance of all 
particle filters when they use 2000 particles; qualitatively similar results were observed for 1000 and 
5000 particles. When the thresholds are chosen so that the proposed filter has the same computational 
complexity as SEPF-EKS, it achieves significantly better tracking performance. Alternatively, for the 
same fixed RMS error performance, the selective processing algorithm reduces the computational load 
by 30 — 40%. Compared to the OOSM-GARP algorithm, a 50% reduction in computational requirements 
leads to only a small increase in estimation error. The results illustrate that we can adjust C ave to control 
the trade-off between the average computational load or power supply consumption and the tracking 
performance. 

VIII. Conclusions 

This paper presents a framework for selective processing of the out-of-sequence measurements. Based on 
this framework we develop a computationally efficient algorithm for delay-tolerant particle filtering that 
has limited memory requirements. By identifying and discarding the uninformative delayed measurements, 
the algorithm reduces the computational requirements. By processing the most informative measurements 
with a re-run particle filter, the algorithm achieves better tracking performance than the storage efficient 
particle filter of lfl~5l . 

In our framework, the threshold to discard uninformative measurements is set by minimizing the one- 
step MSE calculated from the Gaussian approximation of posterior at every filtering time instant. The 
threshold setting could be improved by employing a finite horizon dynamic programming technique to 
take into account the MSE reduction over several forthcoming steps. It is also interesting to explore 
whether the fusion centre can provide feedback to the sensor nodes so that they can locally assess 
measurement informativeness. This would allow sensor nodes to avoid unnecessary energy expenditure 
by discarding uninformative measurements prior to transmission. 
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Fig. 2. RMS vs Complexity from 10 simulations with different C avc - Each simulation shows the average of 1000 MC runs.We 
select three timesteps, t — 10, 20, 30 for filters with 2000 and 5000 particles. The complexity is measured by running time for 
tracking 40s of each filter. The results are run on a Dell laptop with Genuine Intel(R) CPU T2400 1.83GHz, 0.99GB RAM and 
Win-XP OS. 




Appendix A 
Proof of Theorem 1 

We here provide a proof of Theorem 1 . 

We first state a lemma that is employed within the main proof. Denote the spectral radius of a matrix 
by p(0 = maxj |Aj(-)j. The proof of the lemma involves expanding the variational characterization of the 
spectral radius in terms of the blocks of A and applying the Cauchy-Schwarz inequality to each term in 
the expansion. 
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Lemma 1. Let A G R KxK be a block matrix consisting of blocks Ajj G M. Ki * Mi . Then p(A) < 

Proof: Let {xj} and {yj} be two partitionings of vector X G such that X T = [xf ,x^, . . .], 
Xj G M Xl , Vi and X T = [y^, y T , . . .], yj G M Mj , Vj. Write the variational characterization of the spectral 
radius and expand it in terms of blocks: 

IX^AXI 

" (A)= SfW (44) 

= 3S x^x (45) 

<max E ^^ yjl (46) 
x^o X 7 X 

(47) 



Now, for every summand use the Cauchy-Schwarz inequality and take into account the fact that x^Xj/X T X < 
l,Vi: 



|xfA Myj | Ixfx,! 1 / 2 yjA^A,..^' 2 

"^x- ^ ^ [XTX1V2 | X ^x|V2 < 48 > 

| V T A T A ■ -V-IV2 

< max — ^ — —j- . (49) 



Similarly, note the fact that yjyj < X T X,Vj: 



lv T A T A ■ vlV2 lv T A T A- ■vl 1 / 2 

S£ |X^X|V2 ^ 5g |yJ yj -|l/2 (50) 

/ yi A( ; A,,y ; V 2 
= max ^ — ■ . (51) 

This is the variational characterization of p(A| , -Ajj) 1 ' 2 . Substituting it into (l44l completes the proof. 

■ 

Proof of Theorem 1: 

Employing the EKF linear approximation in ( fTOl ) and using the independence of measurement and 
diffusion noises from each other and from the state and independence of and (jh for any m/nor 
s ^ j, we have for m > n and any s,j: 

■ T 

^■y-K = HmFm^Rj^^HJ (52) 
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Note that for m < n R ysK = R^ ayj = Ti s m Rx m x m F^ m H. J n ■ 

Recall that &z k z k is the block-diagonal matrix whose blocks match the diagonal blocks of Hz k z k - 
We now establish finite upper bounds on two spectral radii, p(Bz k z k — R-z k z k ) and p(Rx k z k ~R-z k x k )- 
Throughout the proof we employ the fact that p(C T AC) < p(A)p(C T C) for a square matrix A and an 
arbitrary real matrix C. This follows from the variational characterization of spectral radius: 

, ~ 4 _ |x r C T ACx| / |x T C T ACx| x T C T Cx\ 
p(C AC) = max = max ^ ^— j (53) 

. / |y T Ay| x T CfC X \ 

< max max = = (54) 

x^o \y^o y J y x 7 x / 

< p(A)p(C T C), (55) 
Since all the diagonal blocks of &z k z k — ^z k z k are zero, we have from Lemma [Hand d52l ): 

p(Bz k z k - Rz k z k ) 

^E E[ E E ^( R ^Ry^) 1/2 + E P^yl^YiY^' 2 

m=k—isdS mtk \n=k—£,n^mj£S n , k 3&S m-k ,j^s 

< K£(K£ - 1) maX maX /9(Ry s yjRy:jy s ) 

s,m j^sVn^m m ™ " m 

< - 1) max max max[p(H^F min R x „ x „H^ T H^R x „ x „F^ ^J) 1 / 2 , 

s,m jjtsVnj^m ' 

n(WF Ttv v H s T H J T!v v F T H J? V/ 2 1 

< Kl(Kl- l)maxp(R x „xJmaxp(H^H^ T ) 1 / 2 maxp(H^ T H^) 1 / 2 maxp(F mn F^ n ) 1 / 2 (56) 

n s,m s,m n<m ' ' 

Observing that Rx k z k ~&z k x k = Es,m R ^^ R l fc ^ and recalling <E) we can write: 

< A7max y0 (F fcim F^ m ) /0 (H^ T H^)p(R XmXm ) 2 . (58) 

s,m ' 

Assumptions .A1-.A3 ensure that the bounds in (l56l ) and (I57T ) are finite. 
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We now develop an upper bound for p(R 2 2 — T$zlz k ) 

P( R zlz h ~ B ztz k ) = P^zIzS^za ~ R Zk z k )Bzl Zk ) 

< p{^z\z k )P^z k z k - n ZkZk ) P (B z l Zk ) 

= \^ n (Rz k z k )\^ n (B Zk z k )p(B ZkZk - K ZkZk ). (59) 

Since the eigenvalues of the block-diagonal matrix are the eigenvalues of its blocks we have: \ m i n (Bz k z k ) = 
min s>m A m i n (Ry i =y ? a ). This implies: 

A m in(B2 Z J = 1™ A m i n (H^Rx m I m H^ T + R^ ts ) 

> minA min (R^ ^ ). (60) 

S,m "m^m 

The last inequality holds because (i) for any matrices A and C A m i n (A + C) > A m j n (A) + A m i n (C) and 
(ii) Hf n R,x m x m Hm i s positive semidefinite. 

Similarly, since Hz k z k is a covariance matrix and as such is positive semidefinite we deduce: 

Amm(Rz fe zJ = A min (B 2fc2fc + (Rz k z k - B Zk z k )) (61) 

> max[0, A min (B^ fc2 J + X m in(Rz k z k - B Zk z k )} (62) 

> minA min (R c = iC ^) - p(B Zk z k - ^z k z k ) (63) 

s,m m * 

The last line is valid provided min Sjm A m j n (R^ ^) > p(Bz k z k — ~R-z k z k ), which holds for sufficiently 
large A m i n due to the finite bound derived for the spectral radius in (|59l ). 

We can now derive the following bound on the expression of interest in the theorem, employing the 



relationship tr(-) = 

I tTK Xk z k K z l z Kz k x k - trHx k z k B z l z Kz k x k \ (64) 

= | tr K Xk z k (Rztz k ~ B zlz k )^z k x k I (65) 

< K£p(R Xk z k (Rztz k ~ Bi^jRz^J, (66) 

< K£p(K z l Zk - B z l Zk )p(R Xk z k Rz k x k ), (67) 

< KeX^ n (Kz k z k )^ n (B Zk z k )p(Bz k z k - B, Zk z k )p(Rx k z k Rz k x k ), (68) 

< Kep(R Xk z k 'Rz k x k )p(Bz k z k - ~&z k z k ) (69) 

~ (min Sjm A min (R^^) - p{B Zk z k - Rz k z k )) min Sim A min (R a? sJ ' 
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The finite bounds on the expressions in the numerator lead us to the conclusion that min s m A m i n (R^^) 
oo => | tTTLxAR-zlzRzkX,, ~ tT~R Xk z k Bzl Zk TLz k X k \ -> 0, completing the proof. 

■ 
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